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We construct solutions of the paraxial and Hclmholtz equations which are polynomials in their spatial 
variables. These are derived explicitly using the angular spectrum method and generating functions. Paraxial 
polynomials have the form of homogeneous Hermite and Laguerre polynomials in Cartesian and cylindrical 
coordinates respectively, analogous to heat polynomials for the diffusion equation. Nonparaxial polynomials are 
found by substituting monomials in the propagation variable z with reverse Bessel polynomials. These explicit 
analytic forms give insight into the mathematical structure of paraxially and nonparaxially propagating beams, 
especially in regards to the divergence of nonparaxial analogs to familiar paraxial beams. © 2011 Optical 
Society of America 
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A major aim of optical physics is to find exact math- 
ematical solutions of electromagnetic wave equations, 
which can be identified with observations of the electro- 
magnetic field. We explore here a set of solutions of the 
scalar paraxial and nonparaxial (Hclmholtz) equations 
arising naturally from Taylor series expansions, which 
are polynomials in x and z, x,y and z (Cartesian coor- 
dinates), or R and z (cylindrical coordinates). 

These polynomial beams arise out of Taylor expan- 
sions of optical fields in free space, and describe the op- 
tical amplitude behavior local to the focal plane and 
optical axis, although they diverge at infinity. They 
are therefore a powerful tool in understanding the lo- 
cal structure of optical fields, without concern for global 
properties (which is the domain of Fourier analysis). Ex- 
amples of such analysis by polynomials and Taylor se- 
ries includes the interaction of optical vortices [UH] and 
optical superoscillation [5]. Polynomial beams also illu- 
minate well-known problems in constructing nonparaxial 
counterparts to simple paraxial beams such as Gaussians 
[BHE]- Paraxial polynomials, which are more straightfor- 
ward than nonparaxial, have recently been considered 
in studies of this relationship [HUTU] • We will derive ex- 
plicit forms of polynomial solutions to the paraxial equa- 
tion Vj_V + 2ikd z tfj = and reduced Helmholtz equa- 
tion V 2 ^ + + 2ikd z ip = {e lkz %/j solves the 3D 
Helmholtz equation when ijj solves this latter equation). 
These forms give further insight into the difference be- 
tween paraxial and exact propagation. 

Optical fields with a known analytic form for some 
constant z may readily be represented by power series 
in the spatial variables. The properties of beams based 
on the special functions of mathematical physics, such as 
Gaussian, Bessel and Airy beams, have analytic forms of 
this type. They are usually specified by an initial ampli- 
tude distribution when z = 0, given by a power series in 
Kx or KR, with K a natural inverse length associated 



with the initial field. K is k x for a 2D plane wave, &r for 
a Bessel beam and Wq for a Gaussian beam. Of course, 
power series expansions must converge with respect to 
their variables to make physical sense. 

For instance, a paraxially propagating cylin- 
drical Gaussian beam has the form i))q = 
e -i? 2 /["-o(i+^/^)]/(i + iz/zn), with waist w and 
Rayleigh distance zr = /cuiq/2. In the waist plane z = 0, 
this function is expanded as an exponential power series 
in R/wq- However, for fixed R, the function can also 
be expanded in z around 0, and the appearance of z in 
the denominator of the exponent implies that for R ^ 
there is an essential singularity at z = i%, and the 
power series in z does not converge beyond the Rayleigh 
distance. By comparison, similar series in z of paraxial, 
propagation-invariant beams, such as plane waves and 



Bessel beams = Jo{kR,R)e 



-izkj l /2k 



converge for 



all z (infinite radius of convergence). The nonparaxial 
counterparts of many paraxial beams (i.e. with the same 
initial z = amplitude distribution), such as tJjq, do not 
converge for any z =/= 0. 

We define paraxial polynomials to be solutions of 
the paraxial equation which arc monomials in the ini- 
tial plane, such as x n or R 2n . In Cartesian coordi- 
nates, we write the n th paraxial polynomial p n (x,z), 
with p n (x, 0) = x n . The Fourier transform (in k) of x n , 
is the derivative (5-function i n S^ (k), so we write the 
paraxially propagating polynomial in terms of an angu- 
lar spectrum integral: 
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where we apply, in the second line, integration by parts 
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n times, and in the third line, the substitution n = 

t^/2k/iz. The expression inside square brackets in is 

the generating function of the Hcrmitc polynomials , 
from which it immediately follows that 



Pn(x,z) 
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that is, the paraxial polynomial which is x" when z = 0, 
is a homogeneous polynomial of order n with the coef- 
ficients of Hcrmitc polynomials, each of whose terms is 
made up of powers of x and y/—2iz/k to give the same 
dimension of [Length]™. Since H n {X) only includes odd 
or even powers of X depending on whether n is odd or 
even, only even powers of y/ —2iz/k occur in the parax- 
ial polynomial, which is then genuinely a polynomial in 
x and z (no fractional powers). Generalizing this argu- 
ment, Cartesian paraxial polynomials with initial data 
x n y m are simply the products p n (x, z)p m (y, z). 

Paraxial polynomials in cylindrical coordinates 
Fee ^Pe n (R,z), with vortices of order ±£ on axis 
(with I > 0), are very easily derived by an identi- 
cal argument to that above, with analytic initial data 
R e e ±u ^P £n (R,0) = Rt+^e^*, 
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which is a homogeneous, associated Laguerre polyno- 
mial 11 in R 2 and -2iz/k. Since L e = 1, i? £ e ±l ^ 
by itself solves the paraxial equation (and also the re- 
duced Helmholtz equation). The occurrence of Hcrmitc 
and Laguerre polynomials here is somewhat analogous 
to the appearance of these polynomials with complex 
arguments in elegant Hermite-Gaussian and Laguerre- 
Gaussian beams [12]. 

Paraxial polynomials have been derived by different 
methods and studied before, for example in [3"ll4"ll9]ll"0| . 
The paraxial equation may be thought of as the heat 
equation S7 2 ip — adtij) = with imaginary time t/a = 
—iz/2k, and in this guise the paraxial polynomials are 
known as 'heat polynomials' [13]. Any analytic solution 
of the heat or paraxial equation, with a specified initial 
field distribution at z = expressed as a power scries, 
can be constructed directly by substituting the appropri- 
ate power of x and y or R with the appropriate paraxial 
polynomial, by linearity of the differential equations. 

As described above, the z = amplitude distribution 
of any beam involves a power series in KR or Kx, with K 
an inverse length; usually we think of this as an expan- 
sion in the spatial variable x, or R about 0. However, 
when we replace x n or R 2n with the polynomials ([4|, 
©, the propagated function is formally a power series 
expansion in K about 0. In fact, this is how paraxial 
polynomials occur in free space paraxial beams: func- 
tions such as the Bessel beam ipB or Gaussian beam ipQ, 
expanded about K = 0, must term- by-term satisfy the 
wave equation, since K does not appear in the paraxial 
equation. Each term in K n must therefore be a paraxial 



polynomial of order n; thus, all analytic paraxial beams 
can be expressed as appropriate sums of paraxial poly- 
nomials, and paraxial beams are generating functions for 
paraxial polynomials. The nonparaxial series expansions 
considered in [5] are effectively in K/k, and we consider 
these briefly below. 

We may also think of the propagation of paraxial 
polynomials as the appropriate initial monomial times 
a Gaussian of asymptotically large width wq (limiting 
to a plane wave). In this limit, approaches infinity, 
and the beam close to z = propagates like the poly- 
nomial. The Fourier transform of such a wide Gaussian 
limits to the derivative ^-function in the derivation ([T]) 
above. This approach was used to create fields with de- 
sired geometric properties - knotted optical vortices - 
using polynomial methods [3], which were then embed- 
ded analytically in more physical Gaussian beams with- 
out affecting their local knot topology. 

Nonparaxial polynomials, which solve the reduced 
Helmholtz equation given above, may be constructed by 
modifying the angular spectrum method. The nonparax- 
ial propagator, instead of involving e~ lK z ^ 2k in the an- 
gular spectrum integral (JlJ (or its cylindrical counter- 
part), involves e !fcz (-i+ v /l ~ K2 / fc2 ) ; to which the paraxial 
propagator is the Fresnel approximation. Now, 
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can be thought of as a generating function for monomials 
(y— iz/2k) 2} by expanding in k. Substituting e~ tK z l 2k 
with e *fcz(-i+\/i-K 2 A 2 ) therefore means that the pow- 
ers in z in the paraxial polynomials are replaced by the 
polynomials generated by e i ' £Z (- 1 +V 1 ~ K2 / fe2 ). 

In fact, there are two related polynomial sequences, 
the reverse Bessel polynomials |14[|15) 



ef(Z) = y^Z^ 2 e z K j±1/2 (Z) 



(7) 



where Kj ±l / 2 {Z) is a modified Bessel function of 
half-integer order [TT]. In the literature, 6j~(Z) is 
usually called the 'reverse Bessel polynomial', and 
clearly 0J (Z) = ZO+^Z). The 'Bessel polynomials' 
Z n 9^(l/Z) arc an orthogonal polynomial sequence, but 
&f(Z) themselves are not. 0j (Z) has the explicit form 
9q(Z) = 1, and for j = 1,2,..., 
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The (9 reverse Bessel polynomials arc given by the ex- 



(9) 



ponential generating function [141115 
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Thus, by elementary substitution, the nonparaxial expo- 
nential in the angular spectrum can be expanded 



00 2 j 1 
3=0 ■> 



(-ikz). (10) 



Thus by ([5]) and (|10[) , nonparaxial polynomials can be 
found from the corresponding paraxial polynomial by the 
substitution of all powers of z by 0~ polynomials, 
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This statement is the main result of this Letter. 

Explicitly, the nonparaxial polynomials corresponding 
to Q and © are 
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from the appropriate forms of the Hermite and asso- 
ciated Laguerre polynomials, and the paraxial polyno- 
mials are recovered asymptotically as k — > 00. Non- 
paraxial polynomials can be generated by expanding 
explicit solutions of the reduced Hclmholtz equation 
around K = 0. This was done in [HS] for the plane wave 
"0pw = exp[iKx+iz(\Zk 2 — K 2 — k)] and the Bessel beam 
^ B = e^JtiKR) exp[iz(Vk 2 - K 2 - k)}, although the 
simple forms (fT2"l) . ([TB"]) were not realized. 

The replacement (|11[) is equivalent to the approach 
of Wiinchc [7], in which operators Ti and Ti in z and 
d z were constructed to find the nonparaxial counterpart 
to a paraxial beam. It can be shown that T\(—ikz) n = 
9~ (—ikz), so substitution ([lip is equivalent to the opera- 
tion of Ti . The action of Ti gives a second family of non- 
paraxial polynomials. Since T\ = (1 — ik~ 1 d z )T 2 , these 
polynomials are the same as (|T2"j) and (|T3"|) but with 
replacing 9j . Rather than satisfying the Dirichlet con- 
dition on the beam's amplitude at z = 0, this second 
family satisfies a Neumann-like condition when z = 
(Ref. [7] Eq. (3.30)). Similar operators involving Bessel 
polynomials were considered for the heat and paraxial 
equations in [T6] . 

A general nonparaxial beam ip(x, z),ip(R, z) is there- 
fore represented by a triple power series. On rearranging, 
this gives a series in whose leading-order term is the 
paraxial beam, and later terms are higher-order correc- 
tions [6T7 . The known failure of many such nonparaxial 
power series to converge for z 7^ can be seen explicitly 
using nonparaxial polynomials. Each monomial z 3 in the 
paraxial series expansion becomes a reverse Bessel poly- 
nomial, whose coefficients increase factorially as pow- 
ers of z decrease. For instance, the coefficient for z of a 
cylindrical beam with initial expansion X)^Lo a n{K R) 2n , 
on nonparaxial propagation on-axis (R = 0), from (|13|) 



and © is a series -ikJ2n=o a n (2n)\(-K 2 /k 2 ) n . To con- 
verge, a n must approach quickly enough to defeat this 
factorial growth, which is not the case for Gaussian and 
Airy beams, although it is for Bessel beams. 

The analytic approach suggests another interpretation 
about this divergence of nonparaxial beams: the spectra 
of nonparaxial beams with divergent power series expan- 
sions contain evanescent components which are neither 
forward nor backwards propagating, as their wavevector 
component in z is imaginary. Evanescent plane waves 
solving the reduced Hclmholtz equation have the form 
e iKx+z(VK 2 -k 2 -ik) w ^ x > k. Such waves, however, 

cannot be represented by a power series expansion about 
K = 0, as this is on the opposite side of the branch 
point singularity at K = k. This suggests that nonparax- 
ial power series expansions never formally converge for 
beams whose spectra include evanescent waves, such as 
Gaussian [18] and Airy beams [19] . In such a situation, it 
is appropriate to asymptotically resum the divergent tail 
of the series in K/k, and such an approach has been pro- 
posed jSITTU] . Further work will reveal the relationship be- 
tween nonpar axiality, evanescence, and divergent series, 
and we believe the polynomial solutions constructed here 
will be a useful tool in these investigations. 



References 

1. J. F. Nye, J. Opt. Soc. Am. A 15, 1132-1138 (1998). 

2. M. V. Berry, J. Mod. Opt. 45, 1845-1858 (1998). 

3. M. V. Berry and M.R. Dennis, J. Phys. A 34, 8877- 
(2001). 

4. M. R. Dennis, R. P. King, B. Jack, K. O'Holleran and 
M. J. Padgett, Nature Phys. 6, 118-121 (2010). 

5. Y. Aharonov, F. Colombo, I. Sabadini, D. C. Struppa 
and J. Tollaksen, J. Phys. A 44, 365304 (2011). 

6. M. Lax, W. H. Louisell, W. B. McKnight, Phys. Rev. A 
11, 1365-1370 (1975). 

7. A. Wiinsche, J. Opt. Soc. Am. A 9, 765-774 (1992). 

8. R. Borghi and M. Santarsiero, Opt. Lett. 28, 774-776 
(2003). 

9. A. Torre, J. Opt. 13, 015701 (2011). 

10. R. Borghi, F. Gori, G. Guattari and M. Santarsiero, Opt. 
Lett. 36, 963-965 (2011). 

11. National Institute of Standards and Technology, Dig- 
ital Library of Mathematical Functions (NIST 2010), 
http: //dlmf .nist .gov/ 

12. A. E. Siegman, J. Opt. Soc. Am. 63, 1093-1094 (1973). 

13. D. V. Widder, The heat equation (Academic Press 
1976). 

14. E. Grosswald, Bessel polynomials (Springer 1979). 

15. S. Roman, The umbral calculus (Dover 2005). 

16. L. R. Bragg and J. W. Dettman, Rocky Mountain 
J. Math. 25, 887-917 (1995). 

17. Q. Cao and X. Deng, J. Opt. Soc. Am. A 15, 1144-1148 
(1998). 

18. C. J. R. Sheppard, J. Opt. Soc. Am. A 18, 1579-1587 
(2001). 

19. A. V. Novitsky and D. V. Novitsky, Opt. Lett. 34, 
3430-3432 (2009). 



3 



